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^ ; 1 INTRODUCTION 



ABSTRACT 

We classify orbits of stars that are bound to central black holes in galactic nuclei. 
The stars move under the combined gravitational influences of the black hole and the 
central star cluster. Within the sphere of influence of the black hole, the orbital periods 
of the stars are much shorter than the periods of precession. We average over the orbital 
motion and end up with a simpler problem and an extra integral of motion: the product 
of the black hole mass and the semimajor axis of the orbit. Thus the black hole enforces 
some degree of regularity in its neighborhood. Well within the sphere of influence, 
(i) planar, as well as three dimensional, axisymmetric configurations-both of which 
could be lopsided-are integrable, (ii) fully three dimensional clusters with no spatial 
symmetry whatsover must have semi-regular dynamics with two integrals of motion. 
Similar considerations apply to stellar orbits when the black hole grows adiabatically. 
We introduce a family of planar, non-axisymmetric potential perturbations, and study 
the orbital structure for the harmonic case in some detail. In the centered potentials 
there are essentially two main families of orbits: the familiar loops and lenses, which 
were discussed in Sridhar and Touma (1997, MNRAS, 287, L1-L4). We study the 
effect of lopsidedness, and identify a family of loop orbits, whose orientation reinforces 
the lopsidedness, an encouraging sign for the construction of self-consistent models of 
eccentric, discs around black holes, such as in M31 and NGC 4486B. 
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Following up on a decade of work with ground-based tele- 
scopes (cf. Kormendy 1982, 1987 and Lauer 1983, 1985), 
observations with the Hubble Space Telescope (HST) have 
clarified, and extended our knowledge of the centres of galax- 
ies (cf. Lauer et. al 1995, Gebhardt et al. 1996). The bright- 
ness profiles are generically cuspy, and many correlations 
have been drawn between cusp steepness, isophote shapes, 
luminosities, kinematics, and other nuclear and global prop- 
erties of galaxies. Most, if not all, galaxies might have cen- 
tral black holes (BH), massive enough to power active nuclei 
when there is an adequate supply of fuel (Rees 1990) . About 
a dozen galaxies have been shortlisted as candidates for pos- 
sessing central BHs (cf. Kormendy and Richstone 1995), 
and the case appears strong for our Galaxy, M31, M32, and 
NGC 3115. BH detection based on spectroscopy of stellar 
light requires careful modeling of the distribution and mo- 
tions of stars in the central regions. Set against this back- 
ground of progress, our understanding of the dynamics of 
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nuclear star clusters is quite poor. Gerhard and Binney 
(1985) argued that cusps and BHs will destroy box orbits, 
which are then replaced by chaotic orbits. The chaotic or- 
bits being rounder, it is difficult to construct strongly non- 
axisymmetric equilibria (cf. Merritt 1997, 1998 for a flavour 
of recent work in this area). Sridhar and Touma (1997) 
constructed, cuspy, non-axisymmetric, scale-free discs whose 
potentials are of Stackel form in parabolic coordinates. The 
dynamics in these potentials is fully integrable. A BH could 
be added at the centre without affecting the integrability of 
motion. Discs without a BH support only "banana" orbits. 
The BH stabilizes a box-like family of orbits, the lenses. 
These models, as they questioned the implicit assumption 
that cusps and BHs imply chaos, perked our interest in fur- 
ther investigating dynamics very close to the BH. Here, we 
explore another route to more regular dynamics in galactic 
nuclei. 

The orbits of stars in galactic nuclei are controlled by 
the combined gravitational forces of the central BH, and the 
self gravity of the cluster. Very close to the BH, the orbits are 
nearly Keplerian — stars bound to the BH move on nearly el- 
liptical orbits, whereas very energetic ones follow hyperbolic 
trajectories (we ignore relativistic effects, so our discussion 
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is applicable only to distances beyond several Schwarzschild 
radii from the BH). If the velocity dispersion in the cluster is 
a, we might say that a BH of mass M has a strong influence 
on stellar orbits inside a sphere of radius, ~ GM/a 2 . 
For fiducial values, M ~ 10 s M , a ~ lOOkms -1 , we ob- 
tain rj, ~ 50 pc , a spatial scale that HST with a resolution of 
0."l can easily resolve much farther than the Virgo cluster j^| 
Keplerian ellipses do not precess, and we might expect that 
orbits that are strongly influenced by the BH might support 
asymmetric, even lopsided structures. Such appears to be 
the case for M31 and NGC4486B, both of which have double 
nuclei (Lauer et.al. 1993, 1996); these might be signatures of 
eccentric nuclear stellar discs (Tremaine 1995). A few other 
galaxies with central asymmetric light distributions are dis- 
cussed in Lauer et.al. (1995). As Lauer et.al. (1996) observe, 
" A thorough understanding of the dynamics of eccentric 
disks might allow us to estimate the BH mass directly from 
the disk shape by relating the scale at which the disc sym- 
metry is broken to the hole mass." 

In this paper, we develop a perturbative approach to- 
ward classification of orbits within the BH sphere of influ- 
ence, taking advantage of the super-integrability (i.e. de- 
generacy) of the Kepler problem. Borrowing an averaging 
technique from planetary dynamics, we introduce slow dy- 
namics of precessional motions in § 2, and discuss integra- 
bility of slow dynamics in two, and three dimensional mo- 
tions, for both time independent potentials, as well as the 
case when the mass of the BH grows adiabatically. We in- 
troduce a family of planar potentials that are, in general, 
cuspy, non-axisymmetric and lopsided. The harmonic case 
is particularly suited to simple, analytic treatment, and we 
devote § 3 to an exploration of orbits in non-axisymmetric, 
and lopsided cases; orbit families that might have relevance 
to double nuclei are briefly discussed. Centred anharmonic 
perturbations, with small non-axisymmetry are considered 
in § 4. § 5 is devoted to outlining the limits of applicability of 
the averaging principle, emergence of resonant families, and 
the outbreak of chaos. § 6 provides a summary of our results, 
paying some attention to the applicability of the averaging 
technique to M31 and NGC 4486B. 



2 SLOW DYNAMICS 

Consider a cluster of stars, orbiting around a BH located at 
the origin. The motion of a test star is determined by the 
Hamiltonian, 



H 



GM 



+ $(r) 



(1) 



where r and v are the position vector and velocity of the 
star, M is the mass of the BH, and <E> is the (time indepen- 
dent) mean gravitational potential of the cluster. Within 
the BH sphere of influence, $, by definition, is a small per- 
turbation. Thus the orbits of stars (that are bound to the 

t It should be noted that a itself is often not independent of 
distance from the BH. In this case, a should be taken as the 
dispersion at r = r^, which results in an equation wherein is 
the unknown quantity to be solved for. Perhaps a better estimate 
of the sphere of influence of the BH is the radius at which the 
mass of the cluster equals the mass of the BH. 



BH) may be thought of as Keplerian ellipses that precess 
and deform on time scales that are slow compared to orbital 
times. Hence it is useful to imagine that each star is a slowly 
precessing elliptical ring, with its mass distributed inversely 
proportional to its speed around its orbit — in the planetary 
context, this idea derives from Gauss (see Hagihara (1971) 
for a discussion of Gauss's idea and its applications) . Rauch 
and Tremaine (1996) introduced it to explore relaxation ef- 
fects in star clusters around supermassive BHs. The ellipti- 
cal ring approach emerges naturally in the averaging process 
discussed below. 

When there is a separation of time scales, and one is in- 
terested in the slow evolution, it makes sense to average over 
high frequency variations: thus, time averages of physical 
quantities are equated to their space averages, an idea that 
has its roots in the works of Laplace, Lagrange and Gauss on 
planetary dynamics. Orbit-averaging has also been applied 
to study the evolution of the orbits of comets, perturbed by 
the Galactic tidal field (Heisler and Tremaine 1986). The 
straightforward way to achieve this is to express the prob- 
lem in appropriate action-angle variables, identify the fast 
and (the two) slow angles, and integrate over the fast an- 
gle. Then the conjugate (fast) action is predicted to have 
no secular evolution. Laplace did this for the solar system, 
and concluded that the semi-major axes of the Keplerian 
ellipses of the planets do not evolve secularly. However, this 
is a reasonable approximation for the solar system only over 
times of order 10 4 yr. A major contribution to limiting the 
validity of averaging over orbital motions of the planets are 
mean motion resonances (i.e. "fast- fast" resonances). A nat- 
ural worry is whether a similar phenomenon will afflict stel- 
lar dynamics in galactic nuclei. Let us consider nuclear star 
clusters in the collisionless limit, when there are an infi- 
nite number of stars contributing to a total finite mass for 
the cluster. Then each stellar orbit can be thought of as 
a being smoothly populated with stars. The gravitational 
force exerted by such an orbit on a test star will not dis- 
play time variation on the orbital time scale (which would 
have obtained in the planetary case, when only one planet 
populates the orbit). Thus we expect no mean motion res- 
onances in star clusters around nuclear BHs.Q Resonances 
between fast and slow motions can be important for orbits 
with semi-major axes ~ ru, and will stand in the way of this 
naive averaging procedure. Such resonances are the seeds of 
chaos, which we explore in § 5. 

To apply the averaging principle, it is convenient to 
cast the Kepler problem in terms of the Delaunay action- 
angle variables: we write these as (I , L, L z ;w, g, h), where 
(I,L,L Z ) are action variables, and (w,g,h) are the respec- 
tive conjugate angle variables.^] We list the basic definitions 
below, and refer the reader to text books (c.f. Plummer 1960, 



3 Of course, clusters might have as few as 10 6 stars, so that each 
orbit is not smoothly populated with stars. Hence there can be 
some discreteness noise, which would make the system somewhat 
collisional. 

^ it is unfortunate that the notation employed by planetary dy- 
namicists overlaps so heavily with those commonly used by galac- 
tic dynamicists for other physical quantities. Our notation is non 
standard, but we hope that it minimizes confusion. 
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Goldstein 1985) for the derivation. 



/ = V GMa , where a is the semi-major axis. 
w — mean anomaly; marks time along Keplerian 

orbit — it is at pericentre, and advances 

by 2n in one circuit. 
L — magnitude of the angular momentum; L < I . 
g — angle from the ascending node to the pericentre. 
L z — z-component of angular momentum. 
h — angle from i-axis to ascending node. (2) 

Expressed in the Delaunay variables, the Hamiltonian for 
the Kepler problem (equation [j] with $ set equal to zero) 
assumes the simple form, H^ cp = — (G 2 M 2 /2I 2 ). Thus w 
advances at a uniform rate, 



«0O - (5) 

V at J kep 



9// kop _ (GMf 



dl 



p 



(3) 



whereas the other five Delaunay variables remain constant. 

Let us assume that we have managed to express "!> in 
terms of the new variables. Then the full Hamiltonian, 



Ti 



1 /GM\ 2 

2 \~) 



(4) 



where, with some abuse of notation, we now regard $ as a 
function of all six Delaunay variables. For a general pertur- 
bation, the only conserved quantity is Ti. The angles advance 
at rates given by 

dw _ <9<E> dg _ 9$ dh _ <9<E> 

~dt ~ [ ) + W dt ~ dL> ~dt ~ dT/ ' ' j 

When $ is small, SI is the fastest frequency in the problem; 
thus w varies in time much faster than g and h, and aver- 
aging simply means that we can integrate the Hamiltonian 
of equation (^) over one circuit of w — this is equivalent to 
treating a star as an elliptical ring. The averaged Hamilto- 
nian, 



_i (GMy 

~2 \ I ) 



(6) 



governs the slow dynamics of precessional motions. Since 
TL is independent of w, we recover Laplace's result that I 
is conserved by the slow dynamics. Furthermore, <E> itself is 

another slow integral of motion.|[J] Without further ado, we 
can reach some general conclusions about slow dynamics in 
the region r -C ffc. ■ Below we list our conclusions in order of 
increasing generality: 

• For razor-thin discs, motion is confined to the z = 
plane. Any &(x,y) has the two slow integrals of motion, / 
and $. Hence the slow dynamics for two dimensional poten- 
tials, however non-axisymmetric or lop-sided, is integrable, 
and a straightforward classification of orbits is possible. 

• For axisymmetric $ in three dimensions, L z is an exact 
integral of motion (since $ is independent of h). We now 



II Ti. is, of course, exactly conserved, but this does not give us 
any extra conserved quantities after averaging. 



have three integrals of motion (L z , I and $), and the slow 
dynamics is again fully integrable. 

• In three dimensions, when $ has no symmetry whatso- 
ever, we still have the integrals I and <3>. The slow motion 
can be chaotic, but it is clear that the chaos is limited. 



2.1 Adiabatic growth of the BH 

Averaging is also applicable to time dependent perturba- 
tions, $(r,t), when the time variations are slower than the 
orbital times; as before, I = V GMa is a quasi-invariant. 
In some scenarios of the formation and subsequent growth 
of the BH, its mass increases adiabatically (Peebles 1972, 
Young 1980); conservation of I implies that the semi-major 
axis shrinks in proportion to the growth of M. When 
the growth time is also much longer than the precessional 
timescales, additional (adiabatic) invariants might arise. 
These will, in general, be related to the actions correspond- 
ing to precessional degrees of freedom (L, L z ,g, h). We again 
consider cases with different spatial symmetry. 

• For razor-thin discs, the conserved actions are I and 
<f> Ldg, where the integral is performed at fixed time, over 
one cycle of motion in the L — g plane. 

• For three dimensional, axisymmetric cases, in addition 
to I and § Ldg, we have L z as an exactly conserved quantity. 

• For configurations that have no spatial symmetry, I is 
in general the only conserved quantity, since resonances or 
precessional chaos might destroy the other adiabatic invari- 
ants. 

It is difficult to make further progress without taking into 
account the self-consistent evolution of the cluster potential. 



2.2 Slow planar dynamics 

For razor-thin discs, orbits are restricted to the x — y plane. 
(I, w) determine the semimajor axis and position on the 
orbit, and {L,g) the eccentricity and orientation of the el- 
lipse in the plane. The only difference from the three di- 
mensional case is that we allow L to take both ± values. 
It proves convenient to set w = r\ — esinr;, where rj is 
the eccentric anomaly, and e is the eccentricity given by 
e 2 = [1 — (L/I) 2 ] < 1. If we imagine cartesian coordi- 
nates, (x , y ) , centred at the origin, such that positive x 
is along the major-axis toward the pericentre, and y is 
parallel to the minor-axis, we have x = a(cosr] — e), and 
y = a Vl — e 2 sin r\. The primed coordinates being rotated 
by angle g with respect to the fixed coordinates, (x,y), we 
obtain 

x = a|cosg ^cos?7— \J 1 — P^j — (.sing sin^j- 



y 



a|sin(? ^cosry— \/T— l 2 ^j +£cosg sinr;| 



(7) 



where £ = L/I is a dimensionless angular momentum which 
takes values, — 1 < £ < 1 . 

The frequency of apsidal precession of a test particle 
will depend on both the size (a) and shape (e) of the orbit. 
We wish to associate a characteristic value of this frequency, 
H, with the size (but not the shape) of the orbit. Hence it 
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proves convenient to take /i to be of order the precession fre- 
quency for circular orbits (a somewhat paradoxical notion, 
but the limit is sensible). Then the ratio, 5(1) = fi(I)/Q, is 
a convenient small parameter for slow dynamics. We now 
use a dimensionless time, r = fit, as an appropriate tempo- 
ral measure for slow dynamics. Defining the dimensionless, 
averaged Hamiltonian, 

^ J )=I=/^( 1 -^^H©> (8) 

the equations of motion take the standard form 

id£_ _dH_ . _ dg _dH 

dr~ dg ' 9 ~ dr~ d£ ' 

In fact, much of the qualitative picture of the orbit families 
can be gleaned by studying the contour plot of H in the 
(g,l) plane (for some value of the constant /). 



2.3 A family of planar model potentials 

To develop some analytic understanding of orbital structure, 
we introduce a family of planar potential perturbations that 
allow for both non-axisymmetry and lopsidedness.^] Let us 
define a quadratic form, 



(9) 



2 / j x2 . (V - d2} 2 

q = (x-di) H - 2 , 

and the family of planar potentials, 



$(x,y) = 



sgn(a)$ r Q (q 2 + r\ 
(<&o/2) log (q 2 + r 2 ) , 



. a/2 



a = 0, 



(10) 



(11) 



The potentials depend on the five structural parameters, 
( d, r c , b , c , a ), and a magnitude parameter (<E>o J"o~ a )- The 
potentials are centred at a point that is displaced from the 
BH by d = (di,d 2 ). Potential isocontours have a core of 
radius r c , and constant axis ratio, b. The exponent a is a 
measure of the potential gradient across the contours; we re- 
quire 2 > q > — 1, a range that includes homogeneous cores 
as well as very centrally concentrated mass distributions. 
The potentials are lopsided for non zero d, and are cuspy 
when r c — 0. When both d = and r c = 0, the dynamics 
itself is scale-free, so we may set a — 1 . Lopsided potentials 
are not scale-free, even when r c — 0; this is because the 
lopsidedness sets a scale (d). The magnitude, ($0 '"o"'*) mea_ 
sures the slowness of the dynamics. Specifically, we choose 
units such that fi = (<f> a a /Ir%) = {$ a a /VtlMa r$). As 
discussed above, the small parameter is (/J./S7) oc a < - 1+Q1 - ) ; 
since a > — 1, averaging is a very good approximation for 
small a. 



** In the lopsided cases, the potential exerts a force on the BH, 
so we should imagine that the BH is 'nailed' down. It is possible to 
improve the situation by including cubic and higher order terms 
in the potential, but at the cost of considerably more complexity. 
We decided to forgo this luxury, since our primary goal in this 
paper is to pick out interesting orbits. 



3 PLANAR HARMONIC PERTURBATIONS 

The slow dynamics of harmonic perturbations, for 
which a = 2, can be studied exactly for arbi- 
trary non-axisymmetry, and lopsidedness. Pfenniger and 
de Zeeuw (1989) studied stellar orbits in the combined fields 
of a point mass, and a non-axisymmetric harmonic poten- 
tial. For small energies they report that the orbits are es- 
sentially regular orbits in the Kepler potential of the point 
mass. As we demonstrate in this section, this is not quite 
true; the lenses and loops that emerge from our investiga- 
tions are only locally Keplerian. They are better described as 
precessing and deforming Keplerian ellipses. They point out 
(correctly) the chaos that occurs for orbits with a ~ r^, and 
discuss orbits of much higher energies. These are the orbits 
that will determine the overall shape of a galaxy, and they 
make the interesting point that, if the x and y frequencies 
of the harmonic potential are commensurate, then the high 
energy orbits will tend to stay well away from the centre, 
and hence avoid being strongly perturbed by the BH. We 
do not discuss their work any further, because our interest 
is dynamics within r^. 

The core radius loses dynamical significance in this case, 
so in our explorations of harmonic perturbations, we set r c — 
0: 



d 2 + ^ 

,- 1 52 



(12) 



Substituting for x and y, the expressions given in equa- 
tions ([?]), gives us the perturbation explicitly in terms of 
the Delaunay variables. We then average over w, and obtain 
the following expression for the slow Hamiltonian 



H = 



t 

■di 



4 v 



+ 3y/l-£ 2 (— co S£ ,+ (l 



' j cos 2g 

\ d 2 . 
e ) — sm q 



(13) 



where we have introduced e = (b~ 2 — 1) as a measure of 
non-axisymmetry. The parameter space is 3 dimensional, 
(e, di/a , dz/a ), so that the slow dynamics of even planar, 
harmonic perturbations is quite rich. Below we study a few 
cases. 



3.1 Centred harmonic perturbation: d± = d 2 — 

The simplest case is a centred (di = d 2 = ) perturbation, 
for which we may without loss of generality, set e > . The 
Hamiltonian is 



H 



-|(l + |)^-|(l-^)cos2, 



(14) 



and the equations of motion are, 

I = -f = -§ e( l-*) sin 2g 



OH 

de 



5e „ 

+ y^cos 2g, 



(15) 



tt A term, [(5/2) + (5e/4) + (di/a) 2 + (1 + e)(d 2 /a) 2 ] , has 
been dropped in the slow Hamiltonian, because it is does not 
contribute to slow dynamics. 
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(c) (d) 

Figure 1. Planar, centred, harmonic perturbation: (a) isocon- 
tours of H in the (g,£) plane for e = 0.25, (b) a lens orbit in real 
space corresponding to H = —0.13, (c) a loop orbit in real space 
corresponding to H = —0.47, (d) isocontours of H in the (g, €) 
plane for e = 1 



When e = 0, the full Hamiltonian is axisymmetric and the 
angular momentum is an exactly conserved quantity. The or- 
bits are the rosette-like figures given in standard text books 
on classical mechanics. The averaged description, of course, 
coincides with this picture; H = —3£ 2 /2 , the isocontours 
of which are simply lines of constant £ in the (g,£) plane. 
These loop orbits are now viewed as Keplerian ellipses that 
precess at a rate, g — —3£ . 

When e / 0, the non-axisymmetry gives birth to a new 
family of orbits, the lenses. Figures la shows contour plots of 
H for the case e = 0.25, for which the isocontours of the per- 
turbation have axis ratio b ~ 0.9 . The lenses are parented 
by the short axis orbits, which appear as the stable fixed 
points located at (n/2,0) and (37r/2,0). Figure lb shows a 
lens orbit (for H — —0.13); £ oscillates about zero, whereas g 
librates about the short axis. When g = n/2, it is maximally 
round, and when the pericentre reaches its maximum devia- 
tion from the short axis, the lens elongates to a line ( "radial 
orbit"); its angular momentum now switches sign, and the 
orbit returns to maximal roundness at g — ir/2, but with the 
opposite sign of £ . Unlike lenses, loops have a definite sign 
for £, and g circulates through 27r in one period. The loop in 
Figure lc is a deformed rosette that spends somewhat more 
time near the long axis, as it circulates around the origin. 
The unstable fixed points at (0, 0) and (tt, 0) correspond to 
the long axis orbits. Loops are separated from lenses by the 
separatrix that straddles the unstable fixed points. The sep- 
aratrix orbit, as well as the loops and lenses that are close 
to it spend a lot of time in the vicinity of the long axis. As 




(a) (ft) 

Figure 2. Planar, centred, harmonic perturbation: (a) isocon- 
tours of H for e = 5.0, (b) A long— axis lens orbit, for H = —6.0. 

e increases, lenses occupy an increasingly larger fraction of 
phase space (see Figure Id). 

When g — nix, equations (^) imply that 1 = and 
g = (e — 3)1. Thus, when e = 3, g = 0, tt are lines of 
fixed points in the (g,£) plane. In fact, when e = 3, the 
region occupied by the short-axis lenses has expanded to its 
fullest, squeezing the loops out of existence. This value of t 
corresponds to axis ratio of two (for the isocontours of the 
harmonic perturbation), a case for which the full (unaver- 
aged) Hamiltonian is exactly integrable (see Appendix Al 
for details). For e > 3, loops are completely absent, and all 
of phase space is populated by lenses (see Figure 2a). The 
long-axis orbits now become stable, and parent families of 
long-axis lenses, one of which is shown in Figure 2b. We can 
demonstrate this change in stability of the long-axis orbit, 
by expanding the H of equation (h4) to lowest order about 
the fixed point (0, 0): 

H = ~ 5 r +(e-3)^ + ^g 2 . (16) 

When e < 3, the coefficients of £ 2 and g 2 are of opposite signs 
(the long-axis orbit is unstable) , and have the same sign for 
e > 3 (when the long-axis orbit is stable). We can also see 
that the short-axis orbit remains stable by expanding H 
about (tt/2, 0). With g = tt/2 + 5g, 

H = b -e-{U + 3)^-\e{5gf. (17) 

The coefficients of both £ 2 and (Sg) 2 have the same sign, so 
the short-axis orbit is always stable. 

3.2 Lopsided harmonic perturbation: di ^ 0, 

d 2 = 0, e = 

New orbits emerge when the centre of the perturbation is 
displaced from the BH. The simplest case occurs for e = 0, 
for which the isocontours of the perturbing potential are 
circles. In this case we may, without loss of generality, set 
d 2 = 0, and let di = d ■ The Hamiltonian, 

H = + 3^y/l-P cosg. (18) 

now has a (cos g) term, which distinguishes between orbits 
that are aligned with, from those that are anti-aligned with 
the lopsidedness of the perturbation. We recall that g is the 
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angle from the major (x) axis to the pericentre of the in- 
stantaneous Keplerian ellipse. Thus orbits with g = n arc 
aligned, whereas g = are anti-aligned with the lopsided- 
ness. The only free parameter in H is d/a, which may be 
thought of as a dimensionless measure of lopsidedness. Be- 
low we explore the orbital structure as a function of this 
parameter. 

The panels in Figures 3 show the dependence of the 
isocontours of H on lopsidedness. Changes in the topology 
of isocontours occur when d/a crosses 0.5 (Figure 3a to 3b), 
and again when d/a crosses unity (Figure 3c to 3d). This 
behavior is best understood by following the location and 
stability of the fixed points. The equations of motion are, 

dH 



og a 



dH 



- I 2 sin g 

d £__ 

a ^/T~P 



cos g, 



(19) 



and the fixed points are determined by requiring t — g — . 
There are four fixed points, located at 



9 = 
9 = 

9 = 



(anti-aligned radial orbit); 
(aligned radial orbit); 



0, I = 0, 
7T, I = 0, 

7T, e = ±\/l-^, (aligned loops). (20) 



As Figures 3 show, the anti-aligned, radial orbit is a 
stable fixed point. The stability can be established by ex- 
panding H to lowest order about (0,0): 

H = 3 d -- 3 -(l + d -)f- 3 - d -g> + .... (21) 
a 2 V a J 2 a 

The coefficients of £ 2 and g 2 have the same sign, so the fixed 
point is stable. 

When d < 0.5a (see Figure 3a), the region of anti- 
aligned lenses, shares a border (a "separatrix" which goes 
through the origin of the I — g plane) with rosette-like 
loops. These loops have a definite sense of circulation, anti- 
clockwise or clockwise accordingly as £ is positive or nega- 
tive. As the average value of \l\ increases, we sample this 
family of loops until we reach that rosette which at g = 
has \l\ = 1. This is a critical orbit, beyond which we enter 
the region occupied by the more interesting family of aligned 
loops, a family which we explore further down. As d increases 
to 0.5a, the separatrix, which keeps apart anti-aligned lenses 
from rosettes, grows to its maximum allowable half-width of 
1, squeezing the rosettes out of existence. As can be seen in 
Figure 3b, (for which d/a = 0.6), the rosettes have disap- 
peared. The other fixed points correspond to orbits that are 
aligned with the lopsidedness. The aligned loops have eccen- 
tricity, e = d/a; they begin as circular orbits when d = 0, 
elongate with increasing lopsidedness, reducing to radial or- 
bits when d = a. As discussed above, when d/a < 0.5, the 
aligned loops are separated from the anti-aligned lenses by 
the rosettes. As d/a increases beyond 0.5, and while d/a < 1, 
the family of aligned loops shrinks, and is now separated 
from the anti-aligned lenses by lens-like orbits surrounding 
the aligned radial orbit. All through this variation of d from 
to d = a, the aligned radial orbit is unstable. At d = a, 
the two aligned loops (with either sign of £) merge with the 
aligned radial orbit, which now remains stable for d > a . 
Again, this change of stability of the aligned, radial orbit 
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Figure 3. Isocontours of H for planar, lopsided, harmonic per- 
turbation (a) d = 0.3a, (b) d = 0.6a, (c) d = 0.8a, (d) d = 1.5a. 



may be verified by expanding the Hamiltonian about (7r, 0). 
Setting g — tv + 8g , to lowest order, 

When d < a, the coefficients of £ 2 and {3g) 2 have opposite 
signs (fixed point unstable, but have the same sign for d > a 
(fixed point stable). Hence there are three kinds of fixed 
points: 

1. The anti-aligned radial orbit is stable for all d/a. This 
parents a family of anti-aligned lenses (see Figure 4a). 

2. The aligned radial orbits are stable when d > a, and 
parents a family of aligned lenses (see Figure 4d). 

3. The aligned loop orbits are stable when d < a. These 
parent a family of librating loops, to be discussed in some 
detail below. 

There are no orbits that correspond to anti-aligned 
loops. 



3.3 Possible applications to lopsided galactic 
nuclei 

The aligned loops (which are fixed points of the slow dynam- 
ics) have constant angular momentum; £ is positive (nega- 
tive) for anti-clockwise (clockwise) motion of the particle 
on its Keplerian ellipse. Their eccentricity, e = d/a , which 
means that the centres of the ellipses coincides with the cen- 
tre of the perturbation (the BH is at one of the foci). For a 
fixed value of d, aligned loops of varying sizes (i.e. a), form a 
family of concentric, confocal ellipses. Not only can stars be 
placed on such orbits, but being a non-intersecting family, 
such nested ellipses can also be the streamlines of gas flows. 
Figure 5a shows one such representative set of aligned loops. 
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(a) 



(b) 



Figure 4. Lens orbits for the planar, lopsided, harmonic pertur- 
bation: (a) an anti— aligned lens orbit for a = 1, d = 0.5, (b) an 
aligned lens orbit for a = 1, d = 1.5. 



1))) 




(a) 



(b) 



Figure 5. Loop orbits for the planar, lopsided, harmonic pertur- 
bation: (a) a set of aligned-loops when d = 1. The ten nested 
ellipses have a = 1.2, (1.2) 2 , . . . , (1.2) 10 ~ 6.19, (b) A librating 
loop for d = 0.5, = 1, for which I fluctuates between 0.4 and 
0.97. 



Stable oscillations about these fixed point orbits allows I to 
oscillate, while preserving a definite ± sign. Unlike the cir- 
culating loops of the centred perturbation (see Figure lc), 
g for one of these librating loops oscillates about its mean 
value of 7r (see Figure 5b). 

Tremaine (1995) has suggested that the nucleus of M31 
is an eccentric disc in which stars move on nearly Keplerian 
orbits around a massive BH. In his model, the location of 
the BH coincides with the fainter (brightness) peak, P2, and 
brighter peak, PI, is the apoapsis region of the Keplerian 
ellipses. If we imagine that the lopsided harmonic pertur- 
bation mocks self-gravity, our toy model suggests that the 
orbits such as our aligned loops (and their progeny, the li- 
brating loops) might form the backbone of star clusters such 
as Tremaine's eccentric discs. 



4 CENTRED PERTURBATIONS FOR 
GENERAL a 

When a 7^ 2, it is no longer possible to obtain exact expres- 
sions for the slow Hamiltonian. One possibility is to evaluate 
H numerically, and study the contour plots. Another option 
is to look at limiting cases, which is what we do below. In 
this Section we study centred, nearly circular perturbations: 
di = 0I2 = and 0<e<l. Equation (|ll|), for gives 





(a) 



(b) 



Figure 6. Isocontours of H for e = 0.1 (accurate to first order in 
e): (a) a = 1, (b) a = -0.9. 



us the following scale-free potentials: 



$>(x,y) 



^(^ + (! + ^ 



*/2 



1 I 2 , 

— (r + 



»/2 



(23) 



We can now substitute in equation (|23|), the following ex- 
pressions, 



r — a^l — yl - P cos r/j 
y — dosing ^cos 77 — -\/ 1 — P^j + Icos g sin r\ | , (24) 

for r and y, and average over rj to obtain the averaged Hamil- 
tonian (see Appendix A2 for the derivation) correct to first 
order in e: 



H(e,g;a,e) 



ea ea . 
T - T cos2<?) x 



F ( - 



l+a 



(-01 „2 

+ — t cos 2g x 



Q 
"2 



1-Q 2-Q 



1,1 



2, 1 



(25) 



where F is Gauss' Hypergeometric function. When a = 1 
or 2, the Hypergeometric series for F terminates, and the 
Hamiltonian assumes a simple form. We have already dis- 
cussed the case a — 2 in some detail.^ When a = 1, 

-1/2, 1, 1-f) = 



(26) 



F(0, 1/2,2, 1 
so that 

H(£,g;l,e) = Q + (3 - f) - | (l - I 2 ) cos 2g , (27) 

to first order in e. Figure 6a shows the isocontours of this 
Hamiltonian (a = 1) for e = 0.1 ; as may be seen, these are 
similar to those of Figure la. Keeping e small, we can vary a 
from —1 to 2. Figure 6b shows the isocontours for a — —0.9; 

tt For o = 2, equation ((2^) reduces to equation J14]), which hap- 
pens to be an exact expression, because terms of 0(e 2 ) are iden- 
tically zero for this case. 
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(c) (d) 

Figure 7. Dynamics in the potential of a black hole perturbed 
by a logarithmic potential with b = 0.9: (a) Section on energy 
surface with zero velocity orbit y = and x = OSr^; (b) similar 
section but with x = 2.0r^ ; (c) Periodic orbit at the centre of the 
chain of islands seen in (b); (d) Section on energy surface with 
zero velocity orbit y = and x = 4.0r^. 

again the overall structure is unchanged. Hence we conclude 
that, for small e, the orbital structure is independent of a in 
the range (—1,2), with loops and short-axis lenses broadly 
similar to those in Figures lb and lc. 



5 THE LIMITS OF AVERAGING 

The averaged Hamiltonian is an accurate proxy for the 
full Hamiltonian when the orbital period is much shorter 
than the period of precession of periapse. Such is mostly 
the case for orbits that reside within the sphere of influ- 
ence of the BH. However as one moves outwards, the in- 
fluence of the perturbations relative to the BH's pull in- 
creases, the mismatch in frequencies decreases, thereby in- 
validating the averaging procedure. In practice, the break- 
down occurs at resonances between orbital and precessional 
motion and is responsible for the emergence of minor orbit 
families such as the ones discussed by Gerhard and Bin- 
ney (1985) and Miralda-Escude and Schwarzschild (1989) — 
hereafter GB and MES, respectively. 

In what follows, we display the gradual breakdown of 
averaging as one moves away from a BH of mass M, in the 
scale-free logarithmic potential of MES: 

^ x ,y) = ^ln(x a + ^j, (28) 

where vl is the characteristic velocity of large loops, and b 
is a measure of non-axisymmetry. The dynamics being two 
dimensional, we can explore the phase space structure on a 



Poincare section. Since we are interested in resonances be- 
tween orbital and precessional motions, we choose to strobe 
the dynamics with the orbital period at each apocentre pas- 
sage, recording the angular momentum, L, and the argument 
of the pericentre, g. 

MES studied the scale-free case throroughly. The phase 
space is divided between orbits with a definite sense of cir- 
culation (librating non-aligned loops and rosettes) and those 
with no definite sense of circulation. The latter are broken 
into resonant orbit families: banana, fish, pretzel... etc. The 
resonances are between the orbital period of the star and 
the period of precession of its pericentre. The banana for in- 
stance, performs two orbital revolutions in the time it takes 
its pericentre to precess once. Such resonances are naturally 
missed by the averaged Hamiltonian which assumes that 
the orbital period is much faster than the precessional pe- 
riod (see Touma and Tremaine (1997) for a discussion of the 
averaged Hamiltonian). 

We now place a BH at the origin and study the restruc- 
ting of orbits that occurs. MES did introduce a BH but did 
not show any surface of section. They presented samples of 
high energy orbits, including some regular ones. The BH 
radius of influence is given by, rn = GMjv\. We measure 
distance in units of r^, and look at sections of increasing 
energy. The energy is that of an orbit with zero velocity 
started on the rr-axis, with x = fru, f = 0.4, 2.0, 4.0. 

Within th (Figures 7a), the phase space is occupied by 
two main families of regular orbits: short axis lenses and 
loops. The lenses appear naturally in the integrable cusps 
with BHs presented in Sridhar and Touma (1997). The ap- 
parent integrability of the motion is consistent with the va- 
lidity of averaging at distances where the motion is dom- 
inated by the BH. The phase space of the corresponding 
averaged Hamiltonian was not described in this paper, but 
it hardly differs from the averaged dynamics of centered har- 
monic perturbations discussed above (this should be obvi- 
ous from the results of § 4) . As one moves out in radius past 
Th (Figure 7b), resonances become stronger, and a layer of 
stochastic motion develops near the separatrix. A torus has 
broken up into a chain of islands. The corresponding reso- 
nant periodic orbit, the bowtie, is shown in Figure 7c. Fur- 
ther away from the BH (Figure 7d), the short axis radial or- 
bit becomes unstable, bananas are born, together with other 
resonant orbits families which overlap and lead to large scale 
chaos where once the regular lens orbits lived. 

GB studied the logarithmic potential with core and BH. 
They were mostly interested in energetic orbits far outside 
r% and they find that resonant orbit families and stochas- 
tic orbits dominate the phase space. When commenting on 
dynamics inside r^, they just mention that most orbits are 
loops which, as we now know, is incomplete. In fact, within 
rh, the dynamics is accurately described by the integrable, 
orbit averaged Hamiltonian; the phase space is regular with 
mainly short axis lenses and loops. As one moves outside 
the sphere of influence, lenses and loops break down into 
resonant orbit families, which get stronger, overlap and lead 
to stochastic orbits. The short axis radial orbit goes unsta- 
ble and regular lens orbits disappear, as the region within 
the separatrix is practically engulfed by a single stochastic 
orbit. GB described stochasticity in terms of scattering of 
box orbits by the central cusp and BH. Here, stochasticity 
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emerges with the gradual sacrifice of regular lens orbits to 
strong overlaping resonances as we move away from the BH. 



6 DISCUSSION 

In this paper we have introduced an averaging technique into 
the study of star clusters around massive BHs at the centres 
of galaxies. The dynamics of these clusters is governed by 
the combined gravitational potential of the central BH and 
the cluster itself. Within the sphere of influence (rh) of the 
BH, the cluster potential is a small, yet non negligible per- 
turbation. Orbits that are confined within may be viewed 
as slowly precessing (and deforming) Keplerian ellipses. We 
have shown that these orbits possess a quasi-integral of mo- 
tion, in addition to the energy integral; this nearly conserved 
quantity is I = y/GMa, where M is the mass of the BH, and 
a is the semimajor axis of the orbit (when M is constant, 
the semimajor axis is nearly conserved). The quasi-integral 
arises because the Keplerian orbital frequency is much larger 
than the frequencies of precessional motions. Taking advan- 
tage of this mismatch in frequencies, we averaged over an 
orbital period, to obtain a reduced dynamics of slow, preces- 
sional motions, wherein / appears as a full-fledged integral 
of motion. 

The building of self-consistent stellar models for sys- 
tems with low spatial symmetry (such as the central regions 
of M31 and NGC 4486B, which are known to possess dou- 
ble nuclei) surely requires the existence of more than one 
integral of motion. We now have in hand two integrals of 
motion — the energy and the semimajor axis — that may be 
used to construct dynamical models for the centres of galax- 
ies like M31 and NGC 4486B. Let us verify that the condi- 
tions required for the averaging to apply are met in these 
two cases. 

1. M31 probably has a central BH of mass ~ 6 x 10 7 M . 
For velocity dispersion ~ 160 kms -1 , rh ~ 10 pc (~ 3"). 
The double nuclei are separated by ~ 0."5 arcsec, and much 
of the nuclear disc itself appears to be within r h . Further- 
more, within 1", the mass of the disc and the bulge are 
only 16 %, and 2 % of the mass of the BH respectively (see 
Tremaine 1995). Within r h , the gravitational potentials of 
the disc and the bulge are small, but non negligible, pertur- 
bations to the potential of the BH; indeed, this is a situation 
that is ripe for the application of the averaging technique! 
In § 3.2 and 3.3, we considered orbits in a toy model of a 
lopsided perturbation, and discovered that there are a very 
interesting class of loop orbits that are elongated in the same 
sense as the perturbation. The existence of these aligned 
loop orbits is an encouraging sign that Tremaine's (1995) 
eccentric disc can be realised as a self-consistent model. 

2. NGC 4486B (a low luminosity El companion of M87) also 
has a double nucleus separated by 0."15 (Lauer et.al. 1996). 
There seems to be some evidence for a central BH of mass ~ 
9x 10 8 Mq (Kormendy et. al. 1997). For a velocity dispersion 
~ 175 km s -1 , and an assumed distance to NGC 4486B of 
16Mpc, we estimate that r h ~ l."6, a length scale that 
is not only larger than the separation between the nuclear 
peaks, but also one that is nearly fifteen times larger than 
than the best resolution obtainable with the HST; hence 
the averaging technique may be expected to be useful here 
too. Kormendy et. al. (1997) suggest that the double nucleus 



might arise from a stellar distribution similar to Tremaine's 
(1995) disc. However, as Lauer et.al. (1996) note, there are 
differences from the case of M31, and even the evidence for 
a massive BH is less certain. 

A fundamental result of this paper is that the presence 
of a massive BH enforces a certain degree of integrability 
within its sphere of influence; and this happens because of 
the existence of the secular invariant, I. Even without de- 
tailed knowledge of the (perturbing cluster) potential, we 
were able to note that slow dynamics in (i) razor-thin discs 
is fully integrable because the problem is two dimensional, 
(ii) any axisymmetric potential is fully integrable (note that 
axisymmetry allows for lopsidedness) , and (iii) a potential 
with no spatial symmetry whatsoever still conserves two in- 
tegrals, so that any precessional chaos must be highly limited 
in nature. 

When the potential is time dependent, but with vari- 
ations occurring over times longer than the orbital times, 
averaging is still applicable. We discussed the case of the 
adiabatic growth of the BH. Another case of slow time vari- 
ation comes to mind: this is the 'resonant relaxation' of 
Rauch and Tremaine (1996), a process that relies on two 
body encounters between stars. The time scales are typi- 
cally much longer than orbital times, so averaging applies, 
and the stars in their courses may be treated as elliptical 
rings. Over certain time scales (shorter than the classical 
two body relaxation times), the semimajor axes of all the 
stars are approximately constant, whereas torques between 
the rings exchange angular momentum. As the authors note, 
this can lead to a situation in which some regions of nuclear 
star clusters are relaxed in angular momentum, but not in 
energy. 

We introduced a family of model planar potentials, that 
allow for a range of cusp slopes (a), non-axisymmetry and 
lopsidedness, with a view to classification of orbit families. 
For the harmonic case (a — 2), we showed that centered 
non-axisymmetric perturbations support two main families 
of orbits: short axis lenses and loops. Lenses were introduced 
in Sridhar and Touma (1997) and should be recognized as an 
essential feature of regular motion around BHs. They have 
zero average angular momentum and fill out lens shaped 
regions in configuration space. The loops are the familiar 
rosettes of slightly perturbed Keplerian motion. We explored 
the consequences of lopsidedness and identified a family of 
(resonant) aligned loops whose elongation is in the same 
sense as the lopsidedness. 

When d/2, the case of small non-axisymmetry (e <§C 
1) is analytically tractable. The qualitative nature of the 
orbits (lenses and loops) appears to be insensitive to a. This 
is probably because the force exerted on a star due to the 
perturbing potential — even when it becomes infinitely large 
at the centre — is overwhelmed by the force due to the BH. 
If such is the case for larger non-axisymmetry, then our 
explorations of the harmonic case provide a sketch of the 
possible orbit families in razor-thin discs around massive 
BHs. Whether the qualitative nature of dynamics around 
BHs is indeed independent of the steepness of the cusp in 
density is an unsolved problem, but one that is tractable by 
methods such as our averaging technique. 

Averaging is of course valid as long as we are far away 
from resonances between the orbital period and the peri- 
ods of apsidal or nodal precession. We showed, numerically, 
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that such a condition is indeed satisfied for r < r^, but 
breaks down outside. Orbits which tak e a star far outside 
rw may b e expected to b e chaotic (c.f. ( Gerhard & Binney 
(1985)0 flMerritt (1997)| )) . Central BHs and strong cusps 
might force the main body on an elliptical galaxy to evolve 
toward axisymmetry. However, as we have argued in this pa- 
per, the central regions of galaxies might stubbornly persist 
with their striking variety of non-axisymmetric features. 



in terms of the two functions, A and B: 
A(i;a) = f(-!+£,-2 ,1,1-*° 



B(£;a) 



A 



F 



2 - a 



2 ' 2 



,2, \-t , (B3) 



where we have used formulae 3.666, 8.703 and 9.131 (1) of 
Gradshteyn and Ryzhik (1994). Substituting for A and B in 
equation (p32j) completes the derivation. 
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APPENDIX A: THE CASE e = 3 



We show that even the unaveraged dynamics, in the com- 
bined potential of a BH and a centered harmonic pertur- 
bation, is integrable, when e = 3 (i.e. b — 0.5). Indeed the 
combined potential, 



$ = 



GM 



+ K(x 2 +4y 2 ) 



(Al) 



\/x 2 + y 2 

when expressed in parabolic coordinates (see Sridhar and 
Touma 1997 for details), £ = y + x 2 + y 2 , and r\ — y — 
\J x 2 + y 2 , takes the form 



_ F+(Q F-(r)) 

£ - v v-S, 

where 

*+(€) = K(f-GM) 

F-(T)) = K(—\r]f + GM) . 



(A2) 



(A3) 



This form is that of a separable, Stackel potential in 
parabolic coordinates, which supports integrable motion. 



APPENDIX B: HAMILTONIAN FOR GENERAL 



Let $ denote the potential given in equation (|23|). Substi- 
tute for r and y, the expressions given in equations (p4|), 
and average over 77 in the usual manner. This gives us the 
averaged Hamiltonian, H = ($//i/) = Ho + eHi + 0(e 2 ) , 
where 



Bo 



<j> drj ^1 — \J\ — P cos 77^ 



REFERENCES 

Gebhardt, K., et. al., 1996, AJ, 112, 105 

Gerhard, O. E., & Binney, J., 1985, MNRAS, 216, 467 

Goldstein, H., 1980, Classical Mechanics, 2nd ed. (Reading: 

Addison- Wesley) 
Gradshteyn, I. S., and Ryzhik, I. M., 1994, Table of Integrals, 

Series, and Products, edited by A. Jeffrey (Academic Press, 

Inc, Boston) 

Hagihara, Y., 1971, Celestial Mechanics II, Part I (MIT Press, 

Cambridge, Ma) 
Heisler, J. and Tremaine, S., 1986, Icarus, 65, 13 
Kormendy, J., 1982, in Morphology and Dynamics of Galaxies, 

edited by L. Martinet and M. Mayor (Geneva Observatory, 

Sauverny), p 113 
Kormendy, J., 1987, in Structure and Dynamics of Elliptical 

Galaxies, IAU Symposium No. 127, edited by T. de Zeeuw 

(Reidel, Dordrecht), p 17 
Kormendy, J., & Richstone, D., 1995, ARAA 33, 581 
Kormendy, J., et. al., 1997, ApJ, 482, L139 

Lauer, T.R. , 1983, Ph.D thesis, University of California, Santa 
Cruz 

Lauer, T.R., 1985, ApJ, 292, 104 

Lauer, T.R. et. al. 1993, AJ, 106, 1436 

Lauer, T.R., et. al., 1995, AJ, 110, 2622 

Lauer, T.R., et. al., 1996, ApJ, 471, L79 

Merritt, D., 1997, ApJ, 486, 102 

Mcrritt, D., 1998, Comments on Astrophysics, 19 

Miralda-Escude, J., & Schwarzschild, M., 1989, ApJ 339, 752 

Peebles, P.J.E., 1972, Gen. Rel. Grav., 3, 61 

Pfenniger, D. & de Zeeuw, T., 1989, in Dynamics of Dense Stellar 
Systems, edited by D. Merritt (Cambridge University Press, 
Cambridge), 81 

Plummer, H. C, 1960, An Introductory Treatise on Dynamical 

Astronomy (Dover, New York) 
Rauch, K.P., & Tremaine, S., 1996, New Astronomy, 1, 149 
Rees, M.J., 1990, Science, 247, 817 
Sridhar, S., & Touma, J., 1997, MNRAS 287, L1-L4 
Touma, J. & Tremaine, S., 1997, MNRAS 292, 905 
Tremaine, S., 1995, AJ, 110, 628 
Young, P., 1980, ApJ, 242, 1232 



Hi 



— (p drj ( 1 — yl — £ 2 cos 77 J / 



x I sin g ^cos 77 — \/ 1 — £ 2 ^j + I cos g j . (Bl) 
Some straightforward algebra allows us to express 
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